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A new 8^^0 Phanerozoic database, based on 24,000 low-Mg calcitic fossil shells, yields a prominent 32 Ma 
oscillation with a secondary 175 Ma frequency modulation. The periodicities and phases of these 
oscillations are consistent with parameters postulated for the vertical motion of the solar system across the 
galactic plane, modulated by the radial epicyclic motion. We propose therefore that the galactic motion left 
an imprint on the terrestrial climate record. Based on its vertical motion, the effective average galactic 
density encountered by the solar system is 0.172 + 0.006stat + 0.006sysMopc~^. This suggests the presence 
of a disk dark matter component. 



T 



he Milky Way is a barred spiral galaxy. As a consequence, the solar system revolves around the galaxy and 
carries out a small oscillation in the direction perpendicular to the dense galactic disk, with a period much 
shorter than its orbital period\ For very small amplitudes of oscillation, the disk-crossing period P1/2 

depends solely on the density po in the galactic plane, via P1/2 = (7i;/4Gpo)^^^ = 13.2(po/ [Mopc~^] ) ^^^Ma 
(with G being the gravitational constant, assuming a constant density disk, otherwise the period maybe somewhat 
larger^). Previous determinations of the half-period range between 30 and 42 Ma^"^. The methodologies for 
measuring the density of the galactic disk can be classified into three categories. In the first approach, the known 
components are simply summed up, giving the total baryonic mass density at the galactic plane of about 0.09 to 
0.10 Mopc~^ ^'^ The second approach is based on the vertical kinematics of stars^. Assuming that a given tracer 
population of stars is kinematically relaxed, there is then a relation between the vertical density of the tracer stars 
and their dispersion velocity. This enables derivation of the underlying vertical dependence of the gravitational 
potential, and from it the density, yielding values from 0.07 to 0.26 Mopc~^ ^ l These two approaches are often 
used to estimate the amount of dark matter (DM) in the disk. It was argued that the stellar kinematics are 
inconsistent with a higher end estimate of the central mass density^, implying that DM in the disk may account for 
the deficient mass. Since the Hipparcos data set is presently the most extensive astrometric data available, analyses 
based on it should be considered the most reliable. They appear to converge towards lower values of 0.10 to 
0.13 Mopc~^ suggesting that the unobserved DM in the disk is at most the amount expected from a spherical 
or oblate DM halo, around 0.01 to 0.03 Mopc~^. The DM density is important because its exact value is essential 
for searches of weakly interacting massive particles^. 

The third approach to derive the density of the galactic disk is based on the periodicity of the solar systems 
vertical oscillation (VO). The periodic perturbation of the Gort cloud should increase the population of comets 
crossing the Earth pathway, potentially leaving an imprint in the terrestrial cratering record. Estimates for 
cratering periodicity range from 26 to 36 Ma (e.g., ref. 9), but the sparse statistics are disputed^°. 

Paleo- climatic records could serve as alternative chronometers if a physical mechanism exists to link the solar 
system's vertical motion to the terrestrial climate. Several suggested mechanisms could potentially do so: 

1) perturbation of comets in the Gort cloud that, after disintegration into dust, could potentially cool the climate^ \ 

2) collision with interstellar gas clouds^^, 3) climate modulation via cosmic rays^^"^^ In the latter case, the VG 
would translate into an oscillatory variation in the cosmic ray flux, because the cosmic rays density depends on the 
distance from the galactic plane^^. In addition, the lower interstellar pressure would allow the heliosphere to puff 
up and increase the cosmic ray energy loses as they propagate into the solar system^^. 

Since the cosmic ray mechanism can potentially explain climate variations on longer time scales^^ we will 
assume that the link does indeed operate and with it construct in what follows a model for the average climate. 
However, two points should be noted. First, any VG/climate correlation by itself does not prove that the climate 
link is through any particular mechanism. Second, the cosmic ray climate link is controversial and it should 
therefore not be taken for granted. For example, various criticisms were raised on the validity and implications of 
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Figure 1 | The Gaussian fihered 8^^0 data, separated into four groups. The green Une is the low-latitude, blue the mid-latitude, red the high-latitude, 
and the black line the deep sea subset. The latter three subsets were shifted to minimize the between them and the low-latitude subset (see 
Supplementary Materials for details). Note that the low-latitude data show a warming for the past 15 Ma while the three other subsets exhibit cooling. 
Note also the data gaps around 110 and 210 Ma. The dotted vertical lines denote time intervals used for splicing the different combinations of subsets. 



the long time scale correlations^^, the cloud cover cosmic-ray cor- 
relation on the decade time scale^^ and during Forbush decreases^\ or 
whether the observed link between atmospheric ionization and the 
nucleation of small particles can affect the number density of cloud 
condensation nuclei^^. 

Here we use a new database of b^^O values over 488 Ma, which 
covers almost the entire Phanerozoic, and that transcend the afore- 
mentioned limitations. We will show that it exhibits a very clear 
32 Ma periodicity. 

The data. The original database of some 4500 b^^O measurements 
from low-Mg calcitic fossils^^, a reflection of Phanerozoic climate^^, 
was shown to encompass periodicities of about 145 and 32 Ma^^. 
This database was later updated^^ and once more prior to the 
present study, while keeping its original structure. It now includes 
over 24000 b^^O values. It lists separately data for the deep sea 
(>300 m depth) and surficial waters, the latter separately for the 
low-, mid- and high-latitudes and, where possible, also for the 
subtropical realms. Due to biological evolution, and the resulting 
diversity of fossil phyla, none of the subsets covers the entire 
542 Ma of the Phanerozoic time span. This necessitates splicing of 
the partial records. The entire b^^O and ^^^C database and its 
reduction procedures are described in the Methods below, and 
elaborated in the Supplementary Materials. 

The secular b^^O trend homogenized into 1 Ma bins is reproduced 
in Fig. 1. Visual inspection of this figure, spectral (Supplementary 
Fig. S3) and wavelet analyses (Supplementary Figs. S4, S5) show that 
all subsets (tropical to high-latitude realms) exhibit a similar 30 to 
40 Ma oscillation, confirming its robustness and global validity. This 
permits the splicing procedure over the almost entire Phanerozoic, 
providing corrections for offsets relative to the low-latitude subset, 
the latter covering most of the studied time span, are taken into 
account. We tested 5 possible splicing combinations and they all 
yielded similar outcomes (see Supplementary Fig. S5), with only 
small quantitative differences. All subsequent discussion is therefore 
based on a "master" set that splices the two most complete subsets, 
the mid-latitude one for the <200 Ma interval and the low-latitude 
one for the time interval >200 Ma. The detr ended and high pass 
filtered b^^O data of this "master" (ML200) set are reproduced in 
fig. 2. 

The model. To test whether the apparent b^^O oscillations are related 
to the galactic VO of the solar system, we need to construct a model 



and compare it with the above data. Such a model requires two 
components. One is integrating the galactic orbit, and in 
particular, the VO, while the second component is linking the 
distance from the galactic place to the terrestrial climate. 

For the orbit integration, we assume that there is a negligible drop 
of the density with distance from the galactic plane, implying that the 
vertical oscillation is harmonic. A realistic density fall will typically 
introduce a 10% error^. Namely, the density fitted describes an aver- 
age density experienced by the solar system, not the actual density at 
the galactic plane. We also assume that the density falls exponentially 
with radial distance, such that it can be written as: 



p(r,z)=Poexp(-r/i^). 



(1) 



For an = 3.5 ± 0.5 kpc decay length\ a 650 pc radial epicyclic 
motion of the solar system (RO) would give rise to a 17 ± 3% 
variation in the density (See §S7 for a summary of the kinematic 
parameters). This variation should be higher if any flaring of the disk 
is present. For comparison, the second largest modulation, due to the 
spiral arms, is expected to result in only a 5 to 10% density 
variation^^'^^. 

The orbit can then be integrated to give the VO and a RO, as 
described in Supplementary Materials S3. One of the more interest- 
ing aspects of this orbit is the coupling between the average disk 
density and the vertical oscillation. When the solar system is closer 
to the galactic centre, the density is higher and the vertical oscillation 
is both faster and with a smaller amplitude than when it is further out. 

The main hypothesis in the model is that the average global tem- 
perature depends on the distance from the galactic plane. The lowest 
order expansion of any symmetric relation should be second order, 
that is. 



^^'0 = Ai8+W 



(2) 



The constants will depend on the actual physical link between z and 
the terrestrial temperature, if such a link exists. Again we note that we 
do not assume any particular model to be responsible for this link, 
however, we assume for practical reasons that it is the cosmic ray 
climate link^^ since it provides us with quantitative predictions for 
the actual temperature variations expected. 

We expect the cosmic rays to diffuse out of the galactic plane, 
giving a drop with height. Since the cosmic ray halo is much larger 
than the typical 100 pc amplitude of the VO^^, the lowest order 
expansion should be adequate. Because the oscillations are nearly 
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Figure 2 | The linearly detrended and high pass filtered ML200 6^^0 dataset (in red) for Fourier modes shorter than 49 Ma. This is because we do not 
model processes on longer time scales. The simulated VO motion of the solar system in the galaxy (blue) has a secondary frequency modulation caused by 
the epicyclic motion of the solar system that generates slightly shorter VO periods around 130 Ma and 300 Ma and longer ones in between. 
Because the vertical potential changes adiabatically with the epicyclic motion, the vertical amplitude is larger when the period is longer. The shaded region 
denotes the 95% confidence range for the measured S^^O obtained from the finite number of data points in each bin and the variance in the data. 



harmonic, the T(z) relation associated with eq. 2 would then be close 
to sinusoidal, with an amplitude and frequency that vary with the RO 
oscillations. 

As an example, if the cosmic ray density at the maximum ampli- 
tude Zmax is 30% smaller than at the plane, it would translate into a 2 
± 1°C temperature variation^" (with the plane being colder). This 
would correspond to Bis/zl^^^^0.5±0.25%o without ice-volume 
effects, or more with. 

The free parameters include the present vertical location z and 
velocity v^^qj the average density po and its radial gradient po/R, the 
amplitude of the radial oscillation Ur and time of perigalacticon tp, as 
well as the parameters Ajg and Big appearing in eq. 2. Since many of 
the parameters have independent observations, the best fit para- 
meters in any viable model should be consistent with the observa- 
tions. Note also that there is a degeneracy between z^^^ and Bis, as 
well as between the radial density gradient and the amplitude of the 
radial motion. 

Results 

We optimize these free model parameters by minimizing the resi- 
duals between the predicted model and the detrended "master" set 
(fig. 2). Because the data is noisy, with many local minima, we 
employ a genetic algorithm for minimization. The best fit gives an 
average VO of nearly 32 Ma, with a 175 Ma frequency modulation 
by RO. We use the bootstrap method on the different spliced datasets 
to estimate the statistical and systematic errors, as described in §S4 of 
the Supplementary Materials. For the exact fit parameters see 
Table 1. 

The plot of S^^O as a function of the modeled vertical height z in 
the galactic disk (fig. 3) demonstrates that the average d^^O{z) can be 
adequately described with a quadratic form in z that is assumed in the 
model fit. The figure also shows that these d^^O oscillations are 
present in all latitudinal subsets, again confirming the robustness 
and global validity of the observations. The total S^^O variation is 
about l%o, roughly 4°C, or less if the dynamics of ice caps accounts 
for some of the oscillations in the isotope signal. The impact of the 
secondary modulation of VO by RO, which explains the 30-40 Ma 
range of values around the primary 32 Ma frequency in the wavelet 
analysis (Supplementary Figs. S3-S5), is illustrated in fig. 4. 

Statistical significance of the 32 Ma signal. In order to assess the 
statistical significance of the 32 Ma signal, we adopt the Multitaper 
Method (MTM) often used when studying climate signals. The 
method and the specific toolkit used are described at length in ref 
31. It allows us to estimate the amount of structure at different 
spectral bands and to assess the statistical significance by 



estimating the amount of noise level as a function of frequency 
expected for the null hypothesis. 

Fig. 5 describes the MTM spectrum. Because there is no particular 
preference for any noise model, we conservatively assess the noise 
level as a function of frequency under the assumption of an AR(1) 
modeP\ This modeling allows for more noise to occur naturally at 
low frequencies, ensuring that the significance is not artificially high. 

The spectrum reveals two statistically significant peaks at low 
frequencies, around 150 Ma and around 30 Ma, with a significance 
that is more than 99.9% and 99.99% respectively. If we assume a 
"locally white" model for the null hypothesis instead of an AR(1) 
process, the significance is 99.99% per frequency bin for both peaks. 
(The significance for pure "white noise" will be much higher). 

Since we expect the galactic half period to be 30 to 42 Ma^'^, and 
since the frequency resolution is about 0.001 Ma"\ there are roughly 
10 frequency bins in which the galactic signal can appear. Thus, the 
probability that the given frequency band will randomly yield a signal 
which is as large as that observed in any of the bins, is about 10"^. 

Table 1 reveals that the nominal phases of the geological 32 Ma 
signal and the astronomical data are consistent within 1 . 1 Ma of each 
other (with typical errors of 1.2 to 1.5 Ma on each measurement). 
Thus, there is a about a 15% probability that the two unrelated signals 
will have a phase difference which is as small as the one observed or 
within the combined error. This implies that the probability that 
there will be a random fluctuation giving a galactic like oscillation 
(i.e., within 30 to 42 Ma period) and a phase which is consistent with 
the astronomical data is about 1.5 X lO"'*. 

The probability that this random fluctuation would also have a 
secondary modulation that is consistent both in phase and period 
with the radial epicyclic motion is smaller by at least an order of 
magnitude. Thus the statistical significance of the 32 Ma signal 
and its secondary modulation is of order O(l0~^). 

Discussion 

Given the consistency between the vertical and radial oscillations and 
the paleoclimate data, and the low probability that it could be mim- 
icked by random fluctuations, we conclude with high confidence that 
the terrestrial temperature has a component which is quadratic in the 
distance from the galactic plane. Although this can be naturally 
explained through the cosmic ray climate link^^ the observations 
by themselves do not prove it. 

In addition, it should be noted that although a galactic driver can 
naturally explain a stable —32 Ma cycle, there are terrestrial pro- 
cesses that could drive climate variations on the —32 Ma time scale 
as well. The most prominent is probably mantle convection periodic- 
ally producing plumes that result in large volcanic eruptions/igneous 
provinces. These eruptions will in turn add aerosols and carbon 
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Table 1 | Summary of the paleoclimate derived galactic parameters 

Parameter Paleoclimate based value* Astronomically based value 

VO period 

Density at Epicyclic center: 
the plane Locally: 
VO phase (last plane crossing) 
RO Period 

RO Phase (next perigalacticon) 
Average S^^O range^^ 

*The nominal values are obtained after having shifted the old geological chronology to the recent one^^. The systematic error estimated is the sum of the systematic error obtained by fitting all data 
combinations of (5'^0 subsets and from the updated geological chronology. The statistical error is estimated using a bootstrap method. 

^This value corresponds to the average oscillation period experienced by the solar system, where the epicyclic center is about 0.5 kpc further out than the present location. 

' Different kinematically based density estimates yield a v/ide range of periods, v/ith the recent Hipparcos based determinations giving densities at the lov/er end of this range. See Supplementary Material S6 
for a detailed discussion. 

^Astronomical estimate for the last crossing is based on the solar system being 1 5-30 pc North of the galactic plane, moving northwards at 7 ± 1 km/sec'*'^. 

^ Using Hipparcos based determination of the Oort constants'*'', plus an offset due to the location of the epicyclic center. 

^^The average variation in 5^^0 between maximum warmth and maximum coldness, as obtained in the model fit, corresponds to <4°C. 



3 1 .9 ± 0.63tat ± 0.63ys Ma^ 30-42 Ma * 
0. 1 72 ± 0.006stat ± 0.006sysMopc-^ - 
0.208 ± 0.01 8stat ± 0.01 5sysMopc-^ 0.007-0.25 Mopc" 

3.9 ± 1 .23tat ± O.Bsys Ma 2.8 ± 1 .2 Ma^ 

1 75 ± 93tat ± 1 23y3 Ma 1 72 ± 7 Ma^ 

29 ± 1 73tat ± Bsys Ma 21 ± 7 Ma 

1.03 ± 0.063tat± 0.11 3ys%0 



dioxide to the oceanic- atmospheric system^^ and either cool or warm 
the cUmate. The influence of the Earth mantle connection could also 
drive changes or even reversals in Earth's magnetic field^^'^^, which 
would modulate the atmospheric ionization. The advantage of the 
proposed galactic forcing over a terrestrial driver is that it will pro- 
duce a gradual (sinusoidal) fluctuations as found for most of the d^^O 
record with a relatively steady periodicity, while volcanic forcing will 
more likely produce random abrupt perturbations followed by grad- 
ual relaxations to climate base levels (e.g., "volcanic winters"^^). The 
non-gradual behaviour would also characterize other galactic for- 
cings, such as periodic asteroid/comet bombardments triggered by 
the galactic tides^^. 

As a side note, there is also evidence for a ~ 60 Ma periodicity in 
the literature^^"'*^. For example, Meyers & Peter^^ find a 56 ± 3 Ma 
period in the sedimentation record, which can be explained as an 
orogenic oscillator or as oscillatory dynamics in the mantle convec- 
tion. Thus, these —60 Ma periodicities are probably unrelated to the 
32 Ma cycle discussed here, unless there is a very large north-south 
asymmetry relative to the galactic plane, and the 60 Ma cycle can be 
shown to be in phase with double the 32 Ma cycle discussed here. 

The high accuracy of the estimated VO periodicity also enables 
estimating the density at the galactic plane, and may potentially point 
to an interesting conundrum. Our effective density (for an equivalent 
homogeneous disk) is 0. 172 + 0.0065^^^ + 0.0065^5Mopc^ at the radial 
center of the epicyclic motion, and 0.208 + O.OlSstate + 
0.015sysMopc~^ at the present galactic radius (near perigalacticon, 
see Table 1). These values are inconsistent with the density of 
0.103 + 0.010Mopc~^ obtained by analyzing the dynamics of A & 
F stars in the Hipparcos catalogue^. Some of this inconsistency can be 
explained given that the different numbers do not reflect the same 
density. For example, the different model parameters may have 
changed over the past 500 Ma due to diffusion in parameter space, 
such that the 32 Ma reflects the average witnessed by the solar system 
over the past 500 Ma^^'^'*, and not the exact parameters today. In 
addition, we show in the Supplementary Materials that the 
Hipparcos based determination has a large unaccounted systematic 
uncertainty because the local stellar population is not dynamically 
relaxed. It is large enough to explain the whole inconsistency. It is 
worth to note in this respect that the GAIA star catalogue will be 
published in the near future. Since it will be much more extensive 
than the Hipparcos catalogue, both in terms of number of stars and 
their distance from the solar system, we expect under the above 
model to see convergence between the astronomically derived quant- 
ities and their respective geological counterparts. On the other hand, 
a more extensive geological database could provide additional evid- 
ence in the form of a second "frequency modulation" that corre- 
sponds to the density variations associated with the spiral arm 



passages. These variations are 2-3 times smaller than the RO and 
therefore harder to detect. 

The density at the galactic plane should also be compared to the 
observed baryonic matter and the estimated halo dark matter densi- 
ties, which only contribute 0.09-0.10 Mopc"^ and 0.01-0.03 
Mopc~^ respectively (see ref 4 and references therein). This dis- 
crepancy may imply the presence of a large unaccounted component 
that is not the usual halo dark matter. 

Methods 

Raw data. The first step in the present work was to expand the previously available 
database of S^^O and S^^C measurements for low-magnesium calcite shells of marine 
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Figure 3 | The "master" set 6^^0 data (Fig. 2) plotted against the modeled 
vertical location in the galaxy. Normalization is relative to the maximal Zq 
of the recent maximal excursion around 20 Ma ago. The small orange dots 
are the actual 1 Ma binned data. The thick blue error bars are the averages 
of the combined data binned to 10 equal vertical bins. The red curve is a 
parabolic best fit: S''0/[%o] = 0.42 - 0.09(z/zo) - 0.76(z/zo)l The 
additional points are a coarser binning of the latitudinally separated 
subsets (colours are the same as in Fig. 1). Because of the poor coverage, the 
latitudinal data subsets cannot be high pass filtered. Instead, they have the 
low pass component of the "master" set removed. For legibility, an offset of 
— l%o is applied. The fact that a similar vertical dependence appears in all 4 
unrelated subsets indicates that the phenomenon is a global one. 
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Figure 4 | A raster plot of the detrended 1 Ma binned master set of 8^^0 data. The vertical axis spans the Phanerozoic. The horizontal axis is the 
time folded over a 32 Ma period. For convenience, the horizontal axis is duplicated and shifted by 32 Ma. The blue and red circles (connected by dashed 
lines) are the modeled galactic plane crossings (blue) and the maximal excursions from the plane (red), respectively. An unmodulated 32 Ma signal would 
appear as a vertical line over the entire 490 Ma interval. The apparent modulation due to the radial epicyclic motion (RO) of the solar system is expressed 
as sinusoidal variations of this line. Xs denote bins with insufficient measurements. The disk radii and colour correspond to the detrended and high pass 
filtered S^^O signal, as given by the scale on the right (in %o). 



fossils for the Phanerozoic (last 542 Ma). This present database is expanded by almost 
30% relative to the previous compilation^*^, to about 24,000 measurements (Table SI), 
with brachiopods, foraminifera and belemnites accounting for the bulk of the 
samples. It includes data from more than 20 new references that were published prior 
to July 2011 (see the 2nd sheet of Datafile SI). More information about the 
compilation is available in the Supplementary Materials. The database itself can be 
downloaded as Data File SI. 

Homogenized time series. The first step in the data reduction process was to 
construct temporally homogeneous time series at equidistant 1 Ma intervals. This 



was carried out by Gaussian filtering the raw data of each habitat (tropical to high- 
latitude) separately. The second step was to derive the average S^^O offsets between 
the low altitude and each of the other habitats, by comparing the average S^^O over 
periods with sufficient data in each habitat pair. Because no habitat has a complete 
coverage of the last 490 Ma, the third step was to combine the data sets by splicing 
different habitats together, once the aforementioned offsets are applied to the time 
series. This step produced 5 different sets spanning the last 490 Ma. The different sets 
were all used in the analysis, with the variations giving a handle on possible systematic 
errors. The nominal set was the ML200 which includes the midlatitude data over the 
last 200 Ma, and the low-altitude data before that. It includes only one splice and 
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Figure 5 | The Multitaper Spectrum of the DML75 dataset, using 3 tapers. Without a good model for the noise estimate, the noise level is conservatively 
estimated assuming an AR(1) model, while allowing for more noise at low frequencies than a "white" noise model would. See ref 31 for a description of 
the method. The dominant low frequency peaks are at/~ 0.006 cycle/Ma and/— 0.03 cycle/Ma which are also clearly evident in the wavelet and 
Fourier analyses (see §S2.5 of the supplementary material). Significant higher frequencies are also present. However, they are absent from the wavelet and 
Fourier analyses. It is therefore not clear whether they correspond to real climate variations or to artifacts in the dataset. 
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avoids the recent low-latitude data that show atypical low d^^O values for the last 
15 Ma, if compared to all the habitats. Although the raw data is based on the 2004 
geological chronology, the time scale was adjusted to the new 2012 chronology^^. 
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